function ivp

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0        '
%  where y = (y1, y2, y3 , ..., yn) is an n-vector

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.2,'MR',0.01,'ML',0.05,'P',0.04)

% time interval
tmax=405;

%  initial values
y10=0.1; y20=0;

options = odeset('RelTol',1e-8,'AbsTol',1e-8);
%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],[y10 y20],options);
%[t,y] = ode23s(@rhs,[0 tmax],[y10 y20]);

%
%  2-D plot:  plot solutions as function of t
%
hold on

plot(y(:,1),y(:,2),'r')
plot(y10,y20,'.g','MarkerSize',40)
%plot(t,y(:,2),'--b')
box on
% commands to label each axes
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('y-axis','FontSize',14,'FontWeight','bold')
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14)
hold off


%  define RHS of equations
function dy=rhs(t,y)
b=1;g=1;
dy=zeros(2,1);
dy(1)=y(2)+b*y(1)*y(2)^2;
dy(2)=-y(1)+g*y(1)*y(2)^2;













